/* 3D Navier-Stokes Solver -> Level Set -> Marching Cubes */

#ifndef NAVIER_H
#define NAVIER_H

#include <iostream>
#include <vector>
#include "typen.h"
#include "matrix.h"
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include "velocity.h"

//change the following definition if you have too much objects in your NAV file
typedef unsigned int flagtype;
const flagtype flagtypemsb=UINT_MAX & ~(UINT_MAX >> 1);

/* Class Navier contains all kind of Lists needed for different boundary treatments for
   the velocities, pressure, temperature and chemicals. */
class Navier
{
   protected:

      // Common to other classes - Parameters of the constructor

      int nx;
      int ny;
      int nz;
      double h;

      int size_vectors;

      double **gridUVW;
      double **gridCellUVW;
      int *flag;



      std::vector <int> *List_Border_Face_Air;
      std::vector <int> *List_Border_Face_Fluid;
      std::vector <int> *List_Border_Air;
      std::vector <int> *List_Border_Fluid;

      Scene *S;

      double *gridPhi1;
      double *gridPhi2;

      std::vector <int> **List_Fluid;
      std::vector <int> *List_Fluid_zero;
      int  *List_Fluid_end;
      int  List_Fluid_end_zero;
      int *List_Fluid_sum;
      int *Border_Air_Non_Diagonal_size;
      int *Border_Fluid_Non_Diagonal_size;
      int  Border_Air_size, Border_Air_No_Boundary_size;

      // Computed from the parameters of the constructor
      int dj, dk,hInv;
      double h2;
      double viscosity, froude, tfconv,tfdiff,re,Phi_AIR_FLUID_limit, deltmax;
      double inflow_value[3],g[3],Ext_Force[3];
      int timestepmethod;

      //Parameters changed at every iteration

      double *gridPhi;
      double dt;
      double dt_vector[5];
      double *Accumulated_Output_Flow;



      //Used to define the Initial Scene
      int container_type[6];
      int djx4M4,djx2M2;

      double *Total_inflow;
      double *Total_outflow;
      double *Total_time;


      //Used by Doit routines group
      std::vector <int> List_AirToFluid;
      std::vector <int> List_FluidToAir;
      std::vector <int> List_Cost_Air;
      std::vector <int> List_Cost_Fluid;
      std::vector <int> List_Border_Slip_Outflow_Air;
      std::vector <int> List_Slip_Outflow_Fluid;
      std::vector <int> List_Outflow_Fluid;
      std::vector <int> List_Border_Air_No_Boundary;
      std::vector <int> List_Inflow;

      int  List_Border_Face_size_m1;
      int increment;
      int incrementx2;
      int incrementx4;
      int incrementby2;
      int incrementby4;
      int erase_mask;
      double uvwmax_initial;
      //int size_comp;
      int deltaUVW;

     //Used by Processing_Border_List and New Cell Initializer

     double *output_area[6];

     //created for evaluating the evolution. To be eliminated in the final version
     double uvwmax_NFC_ini;
     double uvwmax;
     double Der_Conv_Max;
     double grid_max_advected[7];
     int    flag_max_advected[7];
     double den_max;
     double sum_coef_max;
     double  dt_max;


      //Used to control the number of iterations
      int num_iter_Navier;
      int iterations, alliterations;//number of iterations for the poisson solver
      int n; //number of steps



      // Used by the advection routines
      double Viscosity;
      double Der_Conv, Laplacian;// variables used for the computation of advection;
      double LL, L, M, R, RR , KL, KR, KRR, KM, KLL, UR, UL, Phi_R, Phi_L;// variables used for the computation of advection;
      bool flag_border_advection; // variables used for the computation of advection;
      int init_FGH,end_FGH;// variables used for the computation of advection;
      double Sum_coef;// variable used for the computation of advection;

      // Used by the Poisson routines
      double *P; // P is the array for pressure
      double  *BiCG_r0, *BiCG_rj, *BiCG_pj, *BiCG_vj, *BiCG_sj, *BiCG_tj, *MatBuf;// Temporary array needed for the BiCGSTAB solver scheme
      double *RHS;// RHS is the right hand side array of the poisson equation
      double norm_res_rjr0_vj_tj;//variavel usada por rotinas do  grupo Poisson.
      double norm_vjr0_tjsj; //variavel usada por rotinas do  grupo Poisson.

      // Used for debugging
      inline void GIJK(int index, int& i, int& j, int& k)
      {
          k =int(index/dk);
          j = (index%dk)/dj;
          i = index%dj;
      }
      //Used for controling precision

      double epsilon;
      double mepsilonxh;

      //Used for updating slip_ouflow_cells
      int d123[3]; // = {1,dj,dk};
      int cell_to_neighbor_delta[6]; //= {-1,1,-dj,dj,-dk,dk};
      int contact_direction[6]; //= {0,0,1,1,2,2};
      int cell_opposite_face_delta[6]; //= {1,0,dj,0,dk,0};
      int cell_opposite_face_delta_0_1[6]; //= {1,0,1,0,1,0};
      int cell_contact_face_delta[6]; //= {0,1,0,dj,0,dk};
      bool  not_contact_table[64][6];
      int index_cell;
      int index_neighbor;
      int neighbor_position;
      int neighbor_code;
      int fluid_neighbors;
      int grid_temp[6];

      //Used for the control of fluid cells taken at each iteration
      int num_iter_min_list_fluid;
      int list_fluid_max;
      int cur_lim_flag_fluid;
      int cur_lim_list_fluid;
      double low_bound_list_fluid;
      double high_bound_list_fluid;
      int  curr_list_fluid_nr;

      double excess_border;


      void Navier_Solver_Initialization();
      void Label_Slip_Outflow_Cells(int index);
      void Init_inflow_and_Label_Slip_Outflow();
      void Call_Update_Slip_Outflow_Cells();
      void Update_One_Neighbor_Slip_Outflow_Cells(bool Outflow_or_Slip);
      void Update_Several_Neighbor_Slip_Outflow_Cells(bool Outflow_or_Slip);
      void Average_several_fluid_neighbors_contribution();
      void Update_Slip_Outflow_Cells(int index,int neighbor);
      void Initialize_Border_Faces_Lists(int index);
      void Processing_Border_Air_List();
      void Ordered_New_Fluid_Cell_Initializer(int curr_index, int curr_antecessors);
      void Diagonal_Border_Cells(std::vector<int>*List,std::vector<int>& List_Cost,std::vector<int>& Target_List, int sign);
      void AdapUVW_Time_step();//velocity correction with pressure gradient due to chorin projection method (see User's Guide)
      void AdapUVW_Time_step_otim();
      void AdapUVW_Border();
      void Call_CompFGH();//calculate temporary velocity fields (in space) due to chorin projection method (see User's Guide)
      void CompFGH(int index);

      double BiCGStab();//BiCGStab poisson solver method
      double BiCGStab_otim();
      double BiCGStab_Border();
      //!calculate temporary velocity fields (in time) due to chorin projection method (see User's Guide
      void CompRHS();// compute right hand side of poisson equation
      //void CompPoisson(/*double starttime*/);//compute poisson equation
      void Poisson();//compute poisson equation
      void Poisson_otim();//compute poisson equation
      void Poisson_Border();
      double CalcRes();
      double CalcRes_otim();
      double CalcRes_Border();
      //  the following functions are used locally
      //interpolate values by third order lagrange interpolation
      //double LagrangeIntp(double LL, double L, double R, double RR, double DD, double D, double P, double PP);
      double MatVecMult(double * A,int index, int RES,  double vjr0_tjsj_factor);//matrix vector multiplication
      double MatVecMult_otim(double * A,int index, int RES,  double vjr0_tjsj_factor);//matrix vector multiplication
       double MatVecMult_Border(double * A,int index, int RES,  double vjr0_tjsj_factor);//matrix vector multiplication
      void  Comp_Der_Conv();


      public:
      void  DoIt();
      void  DoIt_otim();
      Navier(int nx_par,int ny_par, int nz_par, Double h_par,int size_vectors_par,double **gridUVW_par, double **gridCellUVW_par, int *flag_par,
      std::vector <int> *List_Border_Face_Air_par, std::vector <int> *List_Border_Face_Fluid_par,std::vector <int> *List_Border_Air_par, std::vector <int> *List_Border_Fluid_par,
      Scene *S_par, double *gridPhi1_par, double *gridPhi2_par,std::vector <int> **List_Fluid_par,int *List_Fluid_end_par,int *List_Fluid_sum_par,int *Border_Air_Non_Diagonal_size_par,
      int *Border_Fluid_Non_Diagonal_size_par,double *Total_inflow_par, double *Total_outflow_par,double *Total_time_par)
      :nx(nx_par), ny(ny_par), nz(nz_par), h(h_par), size_vectors(size_vectors_par), gridUVW(gridUVW_par), gridCellUVW(gridCellUVW_par), flag(flag_par),
      List_Border_Face_Air(List_Border_Face_Air_par),List_Border_Face_Fluid(List_Border_Face_Fluid_par),List_Border_Air(List_Border_Air_par),List_Border_Fluid(List_Border_Fluid_par),
      S(S_par),gridPhi1(gridPhi1_par),gridPhi2(gridPhi2_par),List_Fluid(List_Fluid_par),List_Fluid_end(List_Fluid_end_par),List_Fluid_sum(List_Fluid_sum_par),
      Border_Air_Non_Diagonal_size(Border_Air_Non_Diagonal_size_par),Border_Fluid_Non_Diagonal_size(Border_Fluid_Non_Diagonal_size_par), hInv(1./h_par),
      Total_inflow(Total_inflow_par),Total_outflow(Total_outflow_par),Total_time(Total_time_par)
     {
        cout<<"Navier Constructor Starting Point \n";
	Accumulated_Output_Flow = new double[size_vectors];
        fill(Accumulated_Output_Flow,Accumulated_Output_Flow+size_vectors,0.0);
	List_Fluid_zero = List_Fluid[0];
        gridPhi = gridPhi1;
	dj = nx + 4;
	dk = (ny + 4)*dj;
	h2 = h*h;
	//viscosity = S->nu;
	viscosity = 0.0;
	froude = S->froude;
        tfconv =  S->tfconv;
        tfdiff = S->tfdiff;
        re = S->re;
	uvwmax_initial = 0;
	for(int jt=0;jt<=2;jt++)
	{
            inflow_value[jt] = S->inflow_value[jt];
	    uvwmax_initial = fmax(uvwmax_initial,fabs(inflow_value[jt]));
	    g[jt]= S->g[jt];
	    Ext_Force[jt]= g[jt]/froude;
	}
	uvwmax = uvwmax_initial;
	Phi_AIR_FLUID_limit = S->Phi_AIR_FLUID_limit;
	timestepmethod = S->timestepmethod;
        deltmax = S->deltmax;
	increment=INCREMENT;
	//increment=256;
        incrementx2 = 2*increment;
	incrementx4 = 4*increment;
	incrementby2 =increment/2;
	incrementby4 =increment/4;
	erase_mask =  ~(6<<3);
	epsilon= S->epsilon;
	mepsilonxh = - epsilon * h;
	*Total_inflow=0.0;
	*Total_outflow=0.0;
	*Total_time=0.0;
	list_fluid_max = S->list_fluid_max;
	low_bound_list_fluid = S->low_bound_list_fluid;
	high_bound_list_fluid = S->high_bound_list_fluid;
	num_iter_min_list_fluid=S->num_iter_min_list_fluid;
        cur_lim_flag_fluid = INFLOW;
        cur_lim_list_fluid = 0;
	cout<<"Ready to initialize List_Fluid_end \n";
	for(int itfor=0; itfor < list_fluid_max; itfor++)
	{
	    List_Fluid_end[itfor] = 0;
	}
	cout<<"Having initialized List_Fluid_end \n";
	List_Fluid_end_zero = 0;
	d123[0]=1;d123[1]=dj;d123[2]=dk;
        int orientation=-1;
        for(int ifor=0;ifor<6;ifor++)
        {
           contact_direction[ifor] = ifor/2; //= {0,0,1,1,2,2}
           cell_opposite_face_delta_0_1[ifor] = 1 - (ifor%2); //= {1,0,1,0,1,0}
           cell_to_neighbor_delta[ifor] = d123[ifor/2]*orientation; //= {-1,1,-dj,dj,-dk,dk}
           cell_opposite_face_delta[ifor] = d123[ifor/2]*cell_opposite_face_delta_0_1[ifor];//= {1,0,dj,0,dk,0}
           cell_contact_face_delta[ifor] = d123[ifor/2] - cell_opposite_face_delta[ifor]; //= {0,1,0,dj,0,dk}
           orientation = -orientation;
        }
	cout<<"Slip_Outflow_Vectors_initialized \n";
	bool  not_contact_table[64][6];

        for(int contact_code=0;contact_code < 64; contact_code ++)
        {
	    int contact_code_divided = contact_code;
            int two_power_face_contact = 32;
	    for(int face_contact=5; face_contact>= 0; face_contact --)
            {
               if((contact_code_divided/two_power_face_contact)==0)
	       {
	           not_contact_table[contact_code][face_contact] = true;
	       }
               else
	       {
	           not_contact_table[contact_code][face_contact] = false;
	       }
	       contact_code_divided = contact_code_divided%two_power_face_contact;
	       two_power_face_contact = two_power_face_contact/2;
	      // cout<<"face_contact,two_power_face_contact ="<<face_contact<<","<< two_power_face_contact<<"\n";
            }
	    // cout<<"contact_code ="<<contact_code<<"\n";
        }
	//cout<<"not_contact_Table constructed \n";
	Navier_Solver_Initialization();
	cout<<"NaSt constructed \n";
     }
      ~Navier();
};

#endif
